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REVISION 2 CHANGES 


This document is a revised version of the previous one dated October 
1971, and supersedes that document. The following is a summary of the major 
technical changes included in this revision: 

1. The filter-weighting matrix extrapolation equations 
have been expanded to provide the (optional) capability 
of including process noise in the filter. 

2. The detailed flow diagrams have been modified so that 
the routine may be re-called to continue an extrapola- 
tion already started on a previous call without re- 
rectification. 

3. The rows and columns of the state covariance matrix 
which pertain to the additionally estimated quantities 
such as landmark locations or instrument biases have 
been re-arranged to be between those rows and columns 
pertaining to the two position and velocity variables. 

This engenders some changes, mostly notational, in 
the filter-weighting matrix extrapolation equations. 
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FOREWORD 


This document is one of a series of candidates for inclusion in a future re- 
vision of MSC-04217, "Space Shuttle Guidance Navigation and Control Design Equa- 
tions". The enclosed has been prepared under NAS9-10268, Task No. 15-A, 

"GN & C Flight Equation Specification Support", and applies to function 1 of the 
Orbital Navigation Module (ON2) and function 1 of the Co-orbiting Vehicle Naviga- 
tion Module (ON3) as defined in MSC-03690 Rev. B, "Space Shuttle Orbiter Guid- 
ance, Navigation and Control Software Functional Requirements - Vertical Flight 
Operations", dated 15 December 1971. 


Gerald M. Levine, Director 
APOLLO Space Guidance Analysis Division 
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NOMENCLATURE 


H d (t) 

b 


c 


nom 


e p 

e t 

f(q) 

G(t) 


— pole 



Hi 


Perturbing acceleration at time t 

Number of additional quantities, such as landmark 
locations or instrument biases, being estimated 

Constant for adjustment of nominal step-size 

Number of columns in the filter weighting sub-matrix W 

Primary vehicle covariance matrix (6x6) 

Target vehicle covariance matrix (6x6) 

Special function of q defined in text 

Gravity gradient matrix 

Unit vector of earth's north polar axis expressed in 
reference coordinates 

Unit vector in the direction of the position vector r 
Three-dimensional identity matrix 

Constant describing dominant term of earth's oblateness 

Special function of £ and 6_ defined in text 

Three-dimensional column vector in the 3x3 process 
noise matrix for either the primary vehicle (i = 3, 4, 5) 
or the target vehicle (i = 9+b, 10+b, 11+b) 
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Q Process noise matrix (3x3). Subscripts P or T 

refer to the process noise matrix associated with 
the primary or target vehicle state 

r 0 Geocentric position vector at time tg 

r(t) Geocentric position vector at time t 

r(t) Magnitude of geocentric position vector 

r (t) Reference conic position vector at time t 

— con 

r„ (t) Magnitude of reference conic position vector 

at time t 

Mean equatorial radius of the earth 
Geocentric position vector at time tp. 

Intermediate values of r 

Switch indicating whether previous extrapolation is 
to be continued without re-rectification 


con 



S cont 


S pert 


Switch indicating the perturbing accelerations to 
be included 


Switch controlling whether process noise is to be 
included in the W-matrix extrapolation 

Switch indicating whether the filter-weighting sub- 
matrix being extrapolated is associated with the 
primary or target vehicle 

Switch controlling whether state or filter-weighting 
matrix integration is being performed (used only 
internally in routine) 
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Xo 


V (t) 

— con 


W, 


W T 


w 


k, i 


w . . 

-k, i 


Initial time point. Also, time of last rectification 

Time to which it is desired to extrapolate (£q. v q ) 
and optionally Wg 

Geocentric velocity vector at time tg 

Geocentric velocity vector at time tp 

Reference conic velocity vector at time t 

Filter-weighting matrix at time tg 

Filter-weighting matrix at time tp 

Elements of the filter-weighting matrix 

Three-dimensional column vectors into which the 
filter- weighting matrix is partitioned 

Independent variable in Kepler routine 

Previous value of x 


y (t) Vector random variable of dimension b representing 

errors in the additionally estimated quantities such 
as landmark locations or instrument biases 

6 (t) Position deviation vector of true position from 

reference conic position at time t 

6' Magnitude of position deviation vector (temporary 

variable used for rectification test) 

6 Maximum value of I 6 | permitted (used as rectifi- 

max — 

cation criterion) 
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At 


Time -step in numerical integration of differential 
equation 


At Maximum permissible time-step size 

At om Nominal integration time-step size 

Time convergence tolerance criterion 

(t) Random variable representing error in estimate of 

position vector at time t 

T7 (t) Random variable representing error in estimate of 

velocity vector at time t 

fi Earth's gravitational parameter 

£(t) Velocity deviation vector of true velocity from refer- 

ence conic velocity at time t 

v' Magnitude of velocity deviation vector (temporary 

variable used for rectification test) 

^ Maximum value of | v | permitted (used as rectifi- 

max 

cation criterion) 

r Time interval since last rectification 

T 1 Previous value of r 

<p Geocentric latitude 
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1 . 


INTRODUCTION 


The Precision State and Filter Weighting Matrix Extrapola- 
tion Routine provides the capability to extrapolate any spacecraft 
geocentric state vector either backwards or forwards in time through 
a force field consisting of the earth's primary central-force gravita- 
tional attraction and a superimposed perturbing acceleration. The 
perturbing acceleration may be either the single dominant term ( J 2 ) 
of the earth's oblateness or a more complete expression involving 
all significant perturbation effects. The Routine also provides the 
capability of extrapolating the filter-weighting matrix along the preci- 
sion trajectory. This matrix, also known as the "W-matrix", is a 
square root form of the error covariance matrix and contains statisti- 
cal information relative to the accuracies of the state vectors and 
certain other optionally estimated quantities. 

On any one call, the routine extrapolates only one state vec- 
tor and only those six rows of the filter- weighting matrix relating to 
this state vector. Two calls are required to extrapolate two separate 
state vectors and a complete filter-weighting matrix pertaining to two 
state vectors. The complete extrapolated filter-weighting matrix is 
obtained by properly adjoining the two separately extrapolated sub- 
matrices of six rows each. 

The routine is merely a coded algorithm for the numerical 
solution of modified forms of the basic differential equations which 
are satisfied by the geocentric state vector of the spacecraft's center 
of mass and by the filter-weighting matrix, namely: 


— £(t) + — r ( t ) 

,,2 3 . . . 

dt r . ( t ) 


= id 


(t) 


and 

~ W(t) = F(t) W(t) + {yQ[ W T (t)] -1 } 

where a^d) is the vector sum of all the desired perturbing accelera- 
tions, F(t) is a matrix containing the gravity gradient matrix and the 
identity matrix in its off-diagonal sub-blocks, and Q is the process 
noise matrix. A simplified form of the term in braces is included 
only during phases when process noise is to be introduced into the 
navigation filter to improve the long-term navigation accuracy. 
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Because of its high accuracy and its capability of extrapola- 
ing the filter-weighting matrix, this routine serves as the computa- 
tional foundation for precise space navigation. It suffers from a 
relatively slow computation speed in comparison with the Conic State 
Extrapolation Routine. 
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2 . 


FUNCTIONAL FLOW DIAGRAM 


The Precision State and Filter Weighting Matrix Extrapola- 
tion Routine performs its functions by integrating modified forms of 
the basic differential equations at a sequence of points separated by 
intervals known as time-steps, which are not necessarily of the same 
size. The routine automatically determines the size to be taken at 
each step. 


As shown in Fig. 1, the state vector and ( optionally) the 
filter-weighting sub-matrix are updated one step at a time along the pre- 
cision trajectory until the specified overall transfer time interval is 
exactly attained. (The size of the last time-step is adjusted as neces- 
sary to make this possible. ) 
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3. 


INPUT AND OUTPUT VARIABLES 


The Precision State and Filter Weighting Matrix Extrapolation Routine has 

the following input and output variables: 

Input Variables 

(r 0< v 0 > Geocentric state vector to be extrapolated. [ If s cQnt f 0, 

(r q , v q ) is last rectified geocentric state vector] 

ig Time associated with (£g,Vg) and W Q . [ If s cont f 0, 

tg is last rectification time] 

tp Time to which it is desired to extrapolate (£ 0 » v Q ) and 

optionally Wq 

s pert Switch indicating the perturbing accelerations to be in- 

cluded. ( s pert = 1 implies J ^ oblateness term only; 

S pert > ^ i m Pli es a more complete perturbing accelera- 
tion model (or models). ) 

d Number of columns in filter- weighting sub-matrix (d = 0, 

6, 7, ... , where 0 indicates no W-matrix extrapolation) 

b Number of additional quantities, such as landmark loca- 

tions or instrument biases, being estimated 

s veh Switch indicating whether the filter-weighting sub-matrix 

being extrapolated is associated with the primary 
(s ve h = 0) or target (s vg ^ = D vehicle 

Wq Filter-weighting sub-matrix to be extrapolated (optional) 

(Wq has dimension 6 x d) 

s Switch indicating whether process noise is to be included 

(s = 1) or not (s = 0) in the W-matrix extrapolation 

q q 

Q Process noise matrix (3x3) associated with the state 

being extrapolated 
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cont 


Switch indicating whether previous extrapolation is to be 

continued (s , = 1) or not (s = 0) without re- 
cont cont 

rectification 


(6,1/) 
(r 
r 


, v ) 
con -con 


(Zp’ — 


Position and velocity deviation vectors 

Conic position and velocity vectors 

Time interval since rectification 

Last value of independent variable in 
Kepler Routine 

Last value of dependent variable in 
Kepler Routine 

Output Variables 
Extrapolated geocentric state vector 


At end of previous 
extrapolation [ used 

onl y if S cont = 1] 


Time associated with (r^, Vp) anb ^F ' 
[ Will equal tp within tolerance of e ] 


Wp Extrapolated filter-weighting sub-matrix of 

dimension 6 x d 


(6_,«i> 


(r , v ) 
— con —con 


^— 0 ’ — 0 ) 


T ' 


Position and velocity deviation vectors 

Conic position and velocity vectors 

Time interval since rectification 

Last rectified position and velocity vectors 

Time of last rectification 

Last value of independent variable in 
Kepler Routine 

Last value of dependent variable in 
Kepler Routine 


For use as input if 

S cont = 1 on a sub " 
sequent extrapola- 
tion 
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4. DESCRIPTION OF EQUATIONS 

4. 1 Precision State Extrapolation Equations 


Since the perturbing acceleration is small compared with the 
central force field, direct numerical integration of the basic differ- 
ential equations of motion of the spacecraft state vector is inefficient. 
Instead, a technique due to Encke is utilized in which only the devia- 
tions of the state from a reference conic orbit are numerically integrated. 
The positions and velocities along the reference conic are obtained 
from the Kepler routine. 

At time t^ the position and velocity vectors, r^ and Vq , define 
an osculating conic orbit. Because of the perturbing accelerations, 
the true position and velocity vectors r (t) and v(t) will deviate as 
time progresses from the conic position and velocity vectors (t) 
and v (t) which have been conically extrapolated from m and v^. 


6 <t) 
^(t) 


r(t) 

v(t) 


^con 

^con 


(t) 

(t) 


be the vector deviations. It can be shown that the position deviation 
6 (t) satisfies the differential equation 


d M 

— 6 (t) + 


dt 


r (t) 
con 


f (q) r (t) + 6 (t) 


= ^(t) 


with the initial conditions 


6_ ( tq ) = 0, v(t Q ) = 0 

where 
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f(q> = q 


( 6_ - 2r ) • 6^ 

q = 

2 

r 


3 + 3q + 

1 + (1 + q) 3 / 2 


an d (t) i- s the total perturbing acceleration. The above 
second order differential equation in the deviation vector <5 ( t ) 
is numerically integrated by a method described in a later sub- 
section. 

The term 

f(q) r(t) +6(t) 

must remain small, i. e. of the same order as a^(t), if the method 
is to be efficient. As the deviation vector 6 (t) grows in magnitude, 
this term will eventually increase in size. When 



I- < t )| > 0.0l|r con (t)| or \u(t)\ > 0. 01 1 ( t ) | 


or when 


I 6 ( t ) I > 6 or I i/( t ) I > i; , 

I — I max I — I max 

a new osculating conic orbit is established based on the latest preci- 
sion position and velocity vectors r(t) and v(t), the deviations 6 (t) 
and r(t) are zeroed, and the numerical integration of 6 (t) and v{t) 
continues. The process of establishing a new conic orbit is called 
rectification. 

The total perturbing acceleration a ^ ( t ) is in general the 
vector sum of all the desired individual perturbing accelerations com- 
prising the total force field, such as those due to the earth's oblate- 
ness, the gravitational attractions of the sun and moon, and the earth's 
atmospheric drag. Since many Shuttle applications will require only 
the perturbing effect of the dominant term J ^ of the earth's oblate- 
ness, the use of only this term has been made a standard option in 
the routine diagrammed in Section 5. However, provision has been 
made for handling a completely general perturbing acceleration. The 
form of this perturbing acceleration will depend primarily upon the 
requirements of the Orbit Navigation function. 
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The explicit expression for the earth's Jg oblateness accel- 
eration alone is: 


id 


JL UL j 
2 2 J 2 


.. 2 r 


(1-5 sin 0 ) _i ^ + 2 sin 0 _ip 0 ^ e 


where 


_i r is the unit position vector in reference coordinates, 


— pole Un ^ vector eart h's north polar 

axis expressed in reference coordinates. 


sin <9 = i • i , 

-r -pole. 


and 


r g is the mean equatorial radius of the earth. 


4. 2 Filter- Weighting (W) Matrix Extrapolation Equations 

The position and velocity vectors which are maintained by the 
spacecraft's computer are only estimates of the actual values of these 
vectors. As part of the navigation technique it is also necessary for 
the computer to maintain statistical information about the position 
and velocity vectors. Furthermore, in particular applications it is 
necessary to include statistical data on various other quantities, such 
as landmark locations during Orbit Navigation and certain instrument 
biases during Co-orbiting Vehicle Navigation. The filter-weighting 
W-matrix is used for all these purposes. 


If e (t) and r) ( t ) are three dimensional vector random 
variables with zero mean which represent the errors in the estimates 
of a spacecraft's position and velocity at time t, then the six-dimen- 
sional state error covariance matrix E(t) at time t is defined by: 


E(t) 


e ( t ) c (t) 


V(t) e ( t ) 


e ( t ) tj( t ) 


V(t) T) ( t) 


where the bar represents the expected value or ensemble average at 
the fixed time t of each element of the matrix over which it appears. 
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If Y(t) is a b-dimensional vector random variable with zero 
mean which represents the errors in the estimates of the b additionally 
estimated quantities such as landmark locations or instrument biases, 
then a ( 6 + b) - dimensional state and other parameter covariance matrix 
is defined by: 



£ (t) £ (t) T 

£ (t) r?(t) T 

e(t) y(t) T 

E(t) = 

2? (t) £(t) T 

h(t) r?(t) T 

2? (t ) y (t) T 


Y(t) e(t) T 

Y(t) T) ( t ) T 

Y(t) y(t) T 


Further, if the statistical properties of the positions and ve- 
locities of two separate spacecraft are to be maintained, a twelve- 
dimensional state covariance matrix is defined by: 


E(t) = 


T 


T 


T 


T 

ip ip 

ip 

ip 

ip 

iT 

ip 

— X 

T 


T 


T 


T 

ip i P 

ip 

ip 

ip 

iT 

iP 

iT 

T 


„ T 


T 


„ T 

-T ip 

iT 

ip 

iT 

iT 

i T 

22, ^ 

T 


T 


T 


T 

iT ip 

ip 

iP 

iT 

iT 

ip 

i T 


where the subscripts P and T refer to the primary and target ve- 
hicles, respectively. 

And finally, if the statistical properties of the b additionally 
estimated quantities are also to be maintained along with the two state 
vectors, a (12 + b) state and other parameter covariance matrix 
is defined by: 
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T 

T 

T 

T 

T 

Ip Ip 

Ip 2p 

I P 2 

Ip It 

f- p 2? 'p 

T 

T 

T 

T 

T 

Sp-p 

22. p 21 p 

RpZ 

22p 1 T 

21 p 22 p 

T 

T 

T 

T 

T 

2 ip 

2 Up 

1 1 

2 Ip 

2 S T 

T 

T 

T 

T 

T 

It -p 

— - T — P 

I t 2 

v rj-* £ rp 

ji. p m. p 

T 

T 

T 

T 

T 

2? T Ip 

Up Up 

H T 2 

22 t — T 

2? T !1 t 


Rather than use one of the above covariance matrices in the 
navigation procedure, it is more convenient to use a matrix W (t) 
having the same dimension d as the covariance matrix E (t) and 
defined by: 

E (t) = W(t) W(t) T 


The matrix W (t) is called the filter-weighting matrix, and is in a 
sense a square root of the covariance matrix. 


Extrapolation of the W (t) matrix in time may be made 
by direct numerical integration of the differential equation which it sat- 
isfies . In the one-spacecraft case, this is: 


_d_ 

dt 



o 

h 1 °(6 x b) 

W(t) = 

G(t) 

° | 


°(b x 6) 

i °(b x b) 


W(t)+ s 1/2 

q 


O 

O 

i o 

1 U (6xb) 
i 

O 

Q 

°(b x 6) 


1 O 

1 (bxb) 


[ W T (t)] 


-1 


(where b = 0, 1, 2, . . . is the number of additionally estimated quantities). 


In the two-spacecraft case, the differential equation is: 




O 

J 3 

1 

o 

G p (t) 

O 

1 

1 


°(b x 6) 


’ 1 
1 
1 

o 

O 

o 

'1 

1 

o 

O 

o 

1 



1 

O 

(6 x b) | 
1 

O 

1 

(b x b) [ 

°(b x 6) 


O 

(6 x b) | 
1 

G s (t) 


o 

o 


o 


W(t) + 
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o 

° i 

°(6 xb) 

I 

1 

O 

O 



o 

Q P ! 

i 

1 

l 

O 

O 


s 1/2 

q ' 

°(b x 6) 

1 

1 

1 

°(b xb) 

1 

1 

i 

°(b x 6) 




O 

o ! 

°(6 xb) 

1 

O 

O 



_o 

° ! 

1 

1 

i 

O 

Q s 



where 1^ is the 3x3 identity matrix, the O's are zero matrices of 
the required dimensions, the G(t ) are the 3x3 conic gravity gradient 
matrices 

G(t) = — — 3r(t)r(t) T - r 2 (t)I 3 
r ( t) L J 

associated with the vehicle under consideration or with the primary (P) 
or target (T) vehicle, and the Q are the likewise associated 3x3 pro- 
cess noise matrices. 

Extrapolation of the W matrix may also be made by the following 
technique, which is somewhat simpler to implement in an on-board computer 
since matrix mainpulations are reduced to more tractable vector manipula- 
tions, and matrix inversion is avoided. 

Let the d x d filter-weighting matrix W = [ w ] be partially 
partitioned into three-dimensional column vectors w , . which bear the 

iC j 1 

subscripts of their first component: 


— o, o—o, r * —o, 5 

l —0, 6 

— 0, 5+b 

1 

1 

— 0, 6+b —0, 7+b 

— 0, d-1 

— 3, 0 — 3, 1 • —3, 5 

j —3, 6 

— 3, 5+b 

1 

1 

1 

— 3, 6+b —3, 7+b 

— 3, d-1 



— 

"T 



W 6,0 W 6, 1 • ’ W 6, 1 

I w 6, 6 

W 6, 5+b 

1 

1 

W 6, 6+b W 6, 7+b 

w 6, d-1 

W 7, 0 W 7, 1* ‘ W 7, 5 

j W 7, 6 

W 7, 5+b 

1 

t 

1 

W 7, 6+b W 7, 7+b 

w 7, d-1 

W 5+b, O' • • w 5+b, 5 

{ W 5+b, 6' ’ 

W 5+b, 5+b 

1 

1 

1 

1 

W 5+b, 6+b W 5+b, 7+b’ * 

' W 5+b, d-1 

— 6+b, O’ ' • — 6+b, 5 

i —6+b, 6' ' 

' — 6+b, 5+b 

1 

1 

i 

— 6+b, 6+b -6+b, 7+b' 

- w 

— 6+b, d-1 

— 9+b, O' ‘ • — 9+b, 5 

1 w 

I —9+b, 6 

• — 9+b, 5+6 

1 

w w 

— 9+b, 6+b — 9+b, 8+b" 

. . w 

— 9+b, d-1 
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and let the 3x3 process noise matrices be partitioned into three-dimensional 
column vectors: 

Qp = [^ 3 ' 5.4 * S 5 1 > Q t = f H 9 +b ’ -3 10 +b’ ^ 1 1 +b ^ 

Furthermore let the inverse of the filter- weighting matrix be approximated by 
the diagonal matrix * (t) whose diagonal elements are the reciprocals of 

diagonal elements of the filter-weighting W-matrix: 

[ W T (t)]" 1 «[ w D (t )] _1 = 


1/w 


0, 0 


1 / W 1, 1 


o 


^d-l, 


d-1 


Then the previous first order differential equations are equivalent to: 


and 


*d t 
with 
w 


2 -0, i 


w 


3, i 
k, i 

j 2 


d t 
d 2 

7^" 

with 


G(t) w Q)i + Sq (l/2w. ( .) fl . 

(i = 3, 4, 5 only) 


dt — 0, i 
constant for k>6 


i = 0, 1 (d-1) 


w 0i G p ( t ) w 0i+ s q (l/2w i(i ) ai 

(i = 3, 4, 5 only) 

— 6+b, i = G T<*> -6+b, i + s q ( + 2 "i, 1’Si 

(i = 9 -t-b, 10+b, 1 1+b only) > i = 0, 1, ... (d- 1) 


w 


3, i 


— 9+b, i 


w 

dt — 0,i 


d t -6+b, i 


w , = constant for 6 2 k < 6+b 

k, 1 


4-7 



(I-P) • • • ‘I ‘0 = I (x-p) ■ 


When written out in full, the above equations are: 

j 2 


T w n ■ = 

d t 2 ~ 0 ’ 1 r •'(t) 


^ j 3 [Kt). w 0 ) .(t) 

+ Sq (1/2 w i £ )£i (i = 3, 4, 5 only) 


r ( t ) - r ( t ) w Q . { t ) 


/ 


with 


w „ . 
— 3, i 


w, . 

k, i 

and 

j 2 


d 

dT 

constant for k a 6 


d t 


2 -0, i 


| 3 [lp (t )‘i 0 ,i(t) 


£p(t) - r p (t) w Qj .(t) 




+ s^ (1/2 w. .)£. (i = 3,4,5 only) 


d 1 2 -8+b.i ' r 5 (t) 


3 [l T (t) -^6+b,i (t) 


r T (t) - r T (t) w 6+b i (t) 


with 


+ s (1/2 w )q (i = 9+b, 10+b, 11+b only) 

M 1, X 1 


w 


3, i 


— 9+b, i 


w 


k, i 


d t -0, i 
d 

d t — 6+b, i 

constant for 6<k<6 + b 


These second-order differential equations may be integrated using the 
same numerical integration technique as is used for the spacecraft 
position vector. The vectors w g . and w 0+b . bear the same relation- 
ship to the spacecraft velocity vector as the vectors w. . and w„ . 

~ ~ U , 1 “ 6^"b, X 

bear to the spacecraft position vector, and w „ . and w . area 

i — 9+b, i 

by-product of the numerical integration of w. . and w_,, . iust as 

— U, l — b+b, l J 

the velocity vector is a by-product of the numerical integration of the 
position vector. 
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4. 3 Numerical Integration Method 

The extrapolation of inertial state vectors and filter weight- 
ing matrices requires the numerical solution of two second -order 
vector differential equations, which are special cases of the general 
form 

,2 

— y(t) = j: (t, y ( t ), z(t)) 
dt 2 


where 


z = — y 
dt “ 


Nystrom's standard fourth-order method is utilized to numerically 
solve this equation. The algorithm for this method is : 


^n+1 = yn + ^ At+ 7<kl + *2 + k 3 )(At)2 

b 

— n+l ‘ i„ + J ( !Si + 2k2 +2 i3 + i4> at 


k, = f (t , y , z ) 
—1 — n* i-n’ -n 


—2 


-3 


*4 


^ (t n + ^ At ' yn + ^ At+ ^ 1 (At) 2 , z^ + ^At) 
*- (t n + ^ At ' Zn + ^ At + Lki(At) 2 , z n + Ik 2 At) 
- (t n + y n + z„ At + l k 3 ( At) 2 , z n + k 3 At) 


where 


y = y (t ), z = z (t ) 
:Z-n i. n " -n — ' n ' 


and 


t , , - t + At 
n+l n 
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As can be seen, the method requires four evaluations of 
£(t, y, z) per integration step At as does the classical fourth-order 
Runge-Kutta method when it is extended to second-order equations. 
However, if £ is independent of z, then Nystrom's method above only 
requires three evaluations per step since = kg. ( Runge-Kutta 's 
method will still require four). 


The integration time step At may be varied from step to 
step. The nominal integration step size is 


At = c r 

nom nom con 


v 2 /jt 


where c is a program constant. (The value c = 0. 3 is 

recommended and implies that about 21 steps will be taken per trajec- 
tory revolution). The actual step-size is however limited to a maxi- 
mum of At max , which is also a program constant. (A value of about 
4000 seconds is suggested. ) Also, in the last step, the actual step 
size is taken to be the interval between the end of the previous step 
and the desired integration endpoint, so that the extrapolated values 
of the state or W-matrix are immediately available. Thus the integra- 
tion step-size At is given by the formula 


At = + minimum (I t„- 1 1, At , At ) 

- '■ F ■’ nom’ max 7 

where tp is the desired integration end-point and t is the time at the 
end of the previous step. The plus sign is used if forward extrapola- 
tion is being performed, while the negative sign is used in the back- 
dating case. 


4-10 



5. 


DETAILED FLOW DIAGRAMS 


This section contains detailed flow diagrams of the Precision State and 
Filter Weighting Matrix Extrapolation Routine. 

Each input and output variable in the routine and subroutine call state- 
ments can be followed by a symbol in brackets. This symbol identifies the 
notation for the corresponding variable in the detailed description and flow dia- 
grams of the called routine. When identical notation is used, the bracketed 
symbol is omitted. 

5. 1 A Note to Coders of the Detailed Flow Diagrams 

The Precision State and Filter Weighting Matrix Extrapolation Routine 
does not require the input of the entire filter weighting matrix. However, for 
coding convenience and conservation of storage, it may be input if desired. If 
only the 6 x d filter weighting submatrix is input, the vehicle switch s ve ^ is 
used only when process noise is included (s^ 0). Its apparent use in Figure 2e 

in order to set the index k is merely so that the notation w . in the flow dia- 
gram will be consistent with the same notation in the description of equations sec 
tion. However, if the entire filter weighting matrix is input, some type of ve- 
hicle switch is necessary even when process noise is not included. The para- 
meters s veh and b could be combined into a single parameter k which is 0 
for the primary vehicle and 6 + b for the target vehicle. For clarity, however, 
they have been kept separate. 
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Figure 2a. Detailed Flow Diagram 


5-2 













Figure 2b. Detailed Flow Diagram 
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(Figure 2e ) 


Figure 2c. Detailed Flow Diagram 
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t = 5. + [j^ + ’ 6 ^ ( k i + k 2 + k 3 ) At At 

v = v + -g <k x + 2k 2 + 2k 3 + k 4 ) At 



Figure 2e. Detailed Flow Diagram 
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SUPPLEMENTARY INFORMATION 


Encke's technique is a classical method in astrodynamics 
and is described in all standard texts, for example Battin (1964). 

The f(q) function used in Encke's technique (and in the lunar-solar 
perturbing acceleration computations) has generally been evaluated 
by a power series expansion; the closed form expression given here 
was derived by Potter, and is described in Battin (1964). 

The oblateness acceleration in terms of a general spherical 
harmonic expansion may be calculated in a variety of ways; three 
different recursive algorithms are given in Gulick (1970). For low 
order expansions, especially those involving mostly zonal terms, an 
explicit formulation is generally superior computation-time-wise, as 
only the non-zero terms enter into the calculation. The general ex- 
pression for the zonal terms is given by Battin (1964), while Zeldin 
and Robertson (1970) give explicit analytic expressions for each of 
the tesseral terms up through fifth order; hence all combinations of 
terms may easily be included in the oblateness acceleration by con- 
sulting the formulations in these references. 

A full discussion of the use of covariance matrices in space 
navigationis given inBattin (1964). Potter(1963) suggested the use of 
the W-matrix and developed several of its properties. It should be noted 
that strictly the gravity gradient matrix G(t) should also include the 
gradient of the perturbing acceleration; however, these terms are so 
small that they may be neglected for our purposes. The use of only 
the conic gravity gradient, however, does not imply the W-matrix is 
being extrapolated conically. (Conic extrapolation of the W-matrix 
can be performed by premultiplying the W-matrix by the conic state 
transition matrix, which can be expressed in closed form). Rather 
the W-matrix is here extrapolated along the precision (perturbed) 
trajectory, as can be seen from the detailed flow diagram of Section 

5. 

The expression for the inclusion of process noise in the 
differential equation satisfied by the filter-weighting matrix is taken 
from Gustafson and Kriegsman (1970), page 7. 
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The Nystrom numerical integration technique was first con- 
ceived by Nystrom ( 1925 ), and is described in all standard texts on 
the numerical integration of ordinary differential equations, such as 
Henrici ( 1962 ). Parametric studies carried out by Robertson (1970) 
on the general fourth-order Runge-Kutta and Nystrom integration 
techniques indicate that the "classic" techniques are the best overall 
techniques for a variety of earth orbiting trajectories in the sense of 
minimizing the terminal position error for all the trajectories, 
although for any one trajectory a special technique can general^ be 
found which decreases the position error after ten steps by one or 
two orders of magnitude for only that trajectory. The classical 
fourth -order Runge-Kutta and Nystrom techniques are approximately 
equally accurate, but the latter possesses the computational advant- 
age of requiring one less perturbing acceleration evaluation per step 
when the perturbing acceleration is independent of the velocity. This 
fact has been taken into account in the detailed flow diagram of Section 
5, in that the extra evaluation is performed only when the perturbing 
acceleration depends explicitly on the velocity. Some past Apollo ex- 
perience has suggested that extra evaluation effect with drag is so 
small as to be negligible; further analysis will confirm or deny this 
for the Space Shuttle. In regard to step-size, the constants and the 
functional form of the nominal and maximum time-step expressions 
have been determined by Marscher (1965). 
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